Especialização em Estatística Aplicada (PPGEAD): Módulo de Estatística Espacial

Introdução à Análise Estatística Espacial

O que é Análise Estatística Espacial ?

São métodos estatísticos que levam em consideração a localização espacial do fenômeno estudado;
Segundo Bailey & Gatrell (1995), “A análise estatística espacial é aplicada quando os dados possuem localização geográfica e quando o arranjo espacial desses dados é considerado relevante para a análise e interpretação dos resultados.”
A primeira questão a ser considerada é: os dados seguem um padrão aleatório ou indicam a presença de agregações bem definidas (clusters) ?

Origem da Estatística Espacial

O uso de dados espaciais na saúde teve um marco histórico com John Snow, que em 1854 mapeou um surto de cólera em Londres, desafiando a teoria miasmática* ao identificar a contaminação da água como causa da doença. Ao mapear as mortes por cólera no Soho, John Snow identificou um foco ao redor da bomba de água da Broad Street, apoiando sua hipótese de transmissão hídrica.

*“Miasmática” refere-se à teoria que defendia que as doenças eram causadas por vapores ou odores nocivos (miasmas) provenientes de matéria orgânica em decomposição, particularmente em áreas úmidas ou com má higiene.

- Mapeamento dos casos de coléra ($\bullet$) e as bombas de água (X) em Londres, 1854.

Dr. John Snow (1813-1858) $\rightarrow$ Considerado pai da Epidemiologia Moderna

A imagem mostra homenagens a John Snow em Soho, Londres: um retrato na fachada de um pub, a réplica da bomba de água da Broad Street, e uma placa que marca a descoberta de que a cólera era transmitida pela água contaminada em 1854.

Objetivos da Estatística Espacial

Investigar padrões espaciais e espaço-temporais, por meio de técnicas como a Análise Exploratória de Dados Espaciais (AEDE) e medidas de correlação espacial, visando identificar estruturas, agrupamentos e dependências nos dados geográficos.
Modelar fenômenos espaciais utilizando modelos estatísticos apropriados, como regressões espaciais (ex: SAR, CAR, GWR) e modelos espaço-temporais, que permitem controlar efeitos de vizinhança (dependência espacial) e heterogeneidade geográfica, com o objetivo de explicar e/ou prever fenômenos influenciados pela localização geográfico.

Dependência Espacial ou Autocorrelação Espacial

Segundo Cressie (1991), embora a suposição de independência entre observações torne a teoria estatística mais tratável, modelos que incorporam dependência estatística costumam ser mais realistas, especialmente em contextos espaciais. Nesse tipo de dado, a dependência entre observações ocorre em múltiplas direções e tende a diminuir conforme aumenta a distância entre os locais amostrados. Em outras palavras, valores próximos no espaço tendem a ser mais semelhantes entre si do que valores distantes, o que caracteriza a autocorrelação espacial.
“Todas as coisas se parecem, porém coisas mais próximas tendem a ser mais semelhantes do que aquelas mais distantes.” (Tobler, 1979). Também conhecida como $1^a$ Lei da Geografia

Tipologia dos Dados Espaciais

Os dados espaciais podem ser classificados segundo diferentes categorias, com base na natureza estocástica de suas observações e na forma como a informação geográfica é representada. Essa tipologia orienta a escolha de métodos estatísticos apropriados para análise.

Segundo Noel Cressie (1993), a estatística espacial pode ser dividida em três grandes áreas:

Dados de Processos Pontuais ou Padrões Pontuais: As observações ocorrem de maneira aleatória no espaço, como casos de uma doença, localização de crimes ou ocorrência de focos de queimadas. O objetivo é entender padrões de agrupamento, dispersão ou aleatoriedade desses pontos.
Dados de Geoestatística: Refere-se as observações que apresentam um atributo mensuravél em localizações contínuas ou irregulares (por exemplo, temperatura, poluição, altitude, teor de argila). Nesses casos, há interesse na dependência espacial entre valores próximos e na interpolação de valores para locais não amostrados, por meio de métodos como krigagem.
Dados de Área: Representam fenômenos agregados por unidades geográficas, como municípios, distritos ou setores censitários. As análises incluem autocorrelação espacial (ex: I de Moran) e modelos de regressão espacial adaptados a dados agregados.

Padrões Pontuais

O principal interesse está no conjunto de coordenadas geográficas representando as localizações exatas de eventos.

Exemplos de aplicação

🏥 Localização de casos de uma doença notificados em uma cidade.
🌳 Distribuição espacial de árvores em um parque urbano.
🐾 Registros de avistamentos de animais silvestres em uma reserva.
🔥 Pontos de ocorrência de focos de incêndio florestal detectados por satélite.

Neste tipo de dado, o evento aleatório de interesse é a posição espacial onde o fenômeno ocorre, e não uma variável medida em si. A análise busca entender se os pontos seguem um padrão aleatório, agrupado (clusters) ou disperso no espaço.

Estimativas de Kernel (ou Mapas de Calor)

As estimativas de densidade kernel são uma técnica amplamente utilizada para representar visualmente a concentração espacial de eventos pontuais, como casos de doenças, crimes ou ocorrências de focos de calor. Essa abordagem transforma um conjunto de pontos em uma superfície contínua de intensidade, evidenciando áreas com maior ou menor densidade (nesse caso a frequencia) de eventos.

Localização da ocorrência de casos de dengue em Belo Horizonte/MG

A técnica foi aplicada para identificar áreas de maior risco de dengue em Belo Horizonte/MG. O mapa resultante revela zonas de alta concentração de casos, o que pode subsidiar ações de controle vetorial e políticas públicas de saúde mais direcionadas.

O objetivo principal é estimar a função de intensidade espacial $\hat{\lambda}_{\tau}(u)$, que descreve a probabilidade de ocorrência de eventos em diferentes locais da região de estudo.

Estimador de intensidade de distribuição de pontos:

\[\hat{\lambda}_{\tau}(u) = \dfrac{1}{\tau^2}\sum k(\dfrac{d(u_i , u)}{\tau}) \text{ , } d(u_i , u) \leq \tau\]

Sabendo que:

$\hat{\lambda}_{\tau}(u)$: estimativa da intensidade do processo pontual no ponto $u$.
$\tau$: parâmetro de suavização (largura de banda), define o raio de influência dos pontos ao redor de $u$.
$u_i$: coordenadas dos pontos observados no plano (ex.: localização dos eventos).
$d(u_i, u)$: distância entre o ponto observado $u_i$ e o ponto de avaliação $u$. Normalmente, a distância euclidiana.
$k(\cdot)$: função kernel que define o peso atribuído a cada ponto com base na distância. Exemplos:

Kernel gaussiano: atribui maior peso aos pontos mais próximos.
Kernel uniforme: todos os pontos dentro do raio recebem o mesmo peso.
Kernel Epanechnikov: pesos decrescem com a distância, chegando a zero na borda.

Condição $d(u_i, u) \leq \tau$: garante que apenas os pontos dentro da vizinhança (de raio $\tau$) contribuam para a estimativa.

Como funciona:

Para cada ponto no espaço, é aplicado um funil de suavização (kernel) que distribui “peso” ao redor do ponto, atribuindo maior peso às áreas próximas e menor peso conforme a distância aumenta.
Ou seja, esse peso decresce com o aumento da distância a partir do centro, de acordo com uma função $k(\cdot)$, e é controlado por um parâmetro chamado largura de banda ($\tau$).
O resultado é uma superfície contínua de densidade, onde regiões com cores mais quentes (vermelho, laranja) indicam maior concentração de eventos.

🌧️ Analogia intuitiva:

O kernel atua como um limpador de para-brisa, pesando mais os eventos próximos e menos os eventos muito distantes. A escolha do tamanho do limpador ($\tau$) e do tipo de movimento (forma da função kernel) define o quanto você vê ao redor do ponto onde está focando.

A figura acima mostra como diferentes larguras de banda ($500 m$, $1.500 m$ e $2.500 m$) alteram o grau de suavização:

Largura pequena ($500 m$): permite detectar agrupamentos muito locais, mas pode introduzir ruídos, ou seja, variações artificiais causadas por flutuações aleatórias nos dados, que não representam padrões reais.
Largura intermediária ($1.500 m$): equilíbrio entre detalhe e suavidade.

Largura maior ($2.500 m$): suaviza padrões locais, ideal para visualizar tendências gerais.

*Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

Geoestatística

Dados usados em geoestatística são atributos contínuos medidos em localizações fixas, na maioria das vezes amostrados no espaço geográfico, e que queremos analisar para entender como um fenômeno varia no espaço.

Estrutura dos dados geoestatísticos

Imagine que você tem um conjunto de pontos no mapa, e para cada ponto você tem um valor observado. O conjunto de dados geralmente tem:

Latitude	Longitude	Atributo Mensurado
-22.90	-43.20	5.4
-22.91	-43.22	6.1
-22.92	-43.18	5.9

Exemplos de aplicação

🌧️ Medição da quantidade de chuva em diferentes locais de uma cidade.
🦟 Contagem de ovos de Aedes aegypti em ovitrampas.
🌫️ Concentração de poluentes no ar em pontos georreferenciados.
🌽 Análise da produtividade agrícola em diferentes talhões de uma fazenda.

Umas das aplicações mais importantes da geoestatística é a interpolação de dados, ou seja, estimação de valores em locais onde não há medição.

Semivariograma: a base da análise espacial

Uma das principais ferramentas da geoestatística é o semivariograma, que mede o quanto dois pontos próximos no espaço se parecem (ou se diferem) com relação ao valor de uma variável.

A ideia central é: Se dois pontos estão muito próximos, é esperado que seus valores sejam parecidos. Já se estiverem distantes, a diferença entre os valores tende a aumentar.

A versão experimental do semivariograma é calculada para diferentes distâncias entre os pontos, usando a seguinte fórmula:

\[\hat{\gamma}(h) = \dfrac{1}{2N(h)}\sum^{N(h)}_{i=1}[z(x_i) - z(x_i + h)]^2 \] Sabendo que:

$\hat{\gamma}(h)$: valor estimado do fenômeno para a distância $h$;
$N(h)$: número de pares de pontos separados por uma distância $h$;
$z(x_i)$: valor observado da variável na localização $x_i$;
$z(x_i + h)$: valor observado da variável em um ponto a uma distância $h$ de $x_i$.
A expressão representa a entre os valores da variável observada em pares de pontos separados pela distância $h$.

📌 Efeito Pepita (C₀): É a variação observada mesmo quando a distância entre os pontos é muito pequena ou zero. Costuma representar erros de medição ou variabilidade microscópica não capturada.

📌 Patamar (C): É o valor em que o variograma se estabiliza, indicando que a partir de certa distância, a variabilidade entre os pontos não aumenta mais. Representa a variância total dos dados.

📌 Alcance (a): É a distância a partir da qual dois pontos deixam de estar correlacionados. Até essa distância, os valores ainda apresentam semelhança espacial, ou seja, após essa distância, tornam-se independentes.

Apesar de útil, essa fórmula não é robusta em todas as situações. Em alguns casos, a variabilidade não é constante ao longo da área estudada — o que chamamos de heterocedasticidade. Nesses cenários, modelos diferentes podem ser utilizados para representar o comportamento do semivariograma. Por exemplo:

Modelo Exponencial
Modelo Esférico
Modelo Gaussiano

Exemplo: Mapa sobre o teor de argila no solo.

A imagem mostra a interpolação do teor de argila em uma área agrícola da Fazenda Canchim (SP), a partir de amostras de solo coletadas em campo. Diversos métodos são comparados: a geoestatística (krigagem), que utiliza a estrutura de dependência espacial dos dados via semivariograma e produz um mapa suavizado e estatisticamente robusto; o método do inverso da distância (IDW), que pondera os pontos mais próximos com maior influência, resultando em maior detalhamento local, mas com risco de exagerar variações; a média simples, que suaviza os dados sem considerar plenamente a variabilidade espacial; e o método do vizinho mais próximo, que gera áreas abruptas ao atribuir a cada região o valor da amostra mais próxima, sem suavização. O exemplo evidencia como a escolha do método de interpolação impacta diretamente a qualidade e continuidade do mapa final.

*Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

Dados de Área

Na análise espacial por áreas, o atributo de interesse costuma ser uma medida agregada (como contagem de casos, taxa de mortalidade, média de renda, etc.) calculada dentro de uma unidade geográfica bem definidas (como bairros, municípios, setores censitários ou regiões administrativas.)
Essas áreas são representadas por polígonos, que podem ter formas irregulares e manter relações espaciais com as áreas vizinhas, seja por fronteiras compartilhadas, conexões físicas (como estradas e rios), ou por semelhanças em características socioeconômicas (como nível de renda ou acesso a serviços).

O objetivo da análise de dados de área é identificar, explicar e interpretar padrões espaciais e tendências que ocorrem entre essas unidades geográficas pré-definidas.

Exemplos de aplicação

🦟 Taxa de incidência de dengue por bairro em uma cidade.
💸 Renda per capita por setor censitário.
🌽 Produtividade agrícola por microrregião.
🏫 Taxa de evasão escolar por município.
🚔 Taxa de criminalidade por distrito policial.
📦 Volume de vendas por zona de entrega.

Mapa Temático

O mapa temático tem como principal objetivo visualizar e analisar a distribuição espacial de um fenômeno específico.
O mapa temático abaixo mostra a frequência de artigos científicos selecionados na revisão integrativa sobre diagnóstico microbiológico de Salmonella spp. na aquicultura entre 2000 e 2020 por países de diferentes continentes. Fonte: Porto et al. 2023 (link).

Matriz de Vizinhança

A matriz de vizinhança $W$ é uma matriz quadrada de dimensão $n \times n$, onde cada elemento $w_{ij}$ representa uma medida de proximidade ou conectividade espacial entre duas regiões $O_i$ e $O_j$.
Essa matriz é fundamental na análise espacial, pois define quais áreas são consideradas vizinhas e como essa vizinhança influencia os fenômenos observados.

A matriz $W$ pode ser binária (0 ou 1) ou ponderada, dependendo do critério adotado.

📌 Critérios comuns para definir vizinhança:

$w_{ij} = 1$, se $O_i$ toca $O_j$ (ou seja, compartilham uma fronteira comum). Esse também é chamado de critério de contiguidade.
$w_{ij} = 1$, se a distância entre os centróides das regiões $O_i$ e $O_j$ for menor que um limite $h$. Esse é o critério baseado em distância geográfica.
$w_{ij} = \frac{l_{ij}}{l_i}$, onde:
- $l_{ij}$: comprimento da fronteira compartilhada entre $O_i$ e $O_j$
- $l_i$: perímetro total da região $O_i$
Esse critério expressa a proporção da fronteira de $O_i$ que está em contato com $O_j$, e é útil para análises mais refinadas de conectividade.

Moran Global, Moran Local e Lisa Map

🧮 Moran Global

Mede a autocorrelação espacial global, ou seja, se valores semelhantes estão agrupados no espaço.

\[I = \frac{ \sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij} (y_i - \bar{y})(y_j - \bar{y}) }{ \sum_{i=1}^{n} (y_i - \bar{y})^2 }\]

$I > 0$: valores similares estão próximos (agrupamento positivo).
$I < 0$: valores diferentes estão próximos (dispersão).
$I \approx 0$: padrão aleatório.

Sabendo que:

$I$: valor do índice de Moran global.
$n$: número total de unidades espaciais (regiões).
$y_i$: valor da variável de interesse na região $i$.
$\bar{y}$: média dos valores de $y$ em todas as regiões.
$w_{ij}$: elemento da matriz de vizinhança que representa a relação espacial entre as regiões $i$ e $j$.
- $w_{ij} = 1$ se $i$ e $j$ são vizinhos, 0 caso contrário (ou ponderado).
O numerador mede a covariância espacial ponderada, e o denominador é a variância total de $y$.

🧮 Moran Local (LISA)

Avalia a autocorrelação espacial local, permitindo identificar regiões com agrupamentos (clusters) ou comportamentos atípicos.

\[I^{(i)} = \frac{n}{\sum_{j=1}^{n} (z_j - \bar{z})^2} \sum_{j=1}^{n} w_{ij} (z_i - \bar{z})(z_j - \bar{z})\]

Sabendo que:

$I^{(i)}$: índice de Moran local para a região $i$.
$n$: número total de regiões.
$z_i$: valor padronizado (ou original, dependendo da convenção) da variável na região $i$.
$\bar{z}$: média da variável $z$.
$w_{ij}$: peso espacial entre as regiões $i$ e $j$.
A soma considera a influência dos vizinhos $j$ sobre o ponto $i$, ponderada pela matriz de vizinhança.

🗺️ Resumindo

Moran Global indica se há padrão geral de autocorrelação espacial (positivo, negativo ou aleatório).
Moran Local revela onde esses padrões ocorrem: clusters de alto/alto (hotspots), baixo/baixo (coldspots), alto/baixo (outliers), etc.

🗺️ LISA Map

O LISA Map é um mapa temático que representa graficamente os resultados do Índice de Moran Local. Ele identifica:

Clustering local (agrupamento de valores semelhantes)
Outliers espaciais (valores discrepantes em relação aos vizinhos)
Regiões não significativas (sem padrão espacial detectável)

Tipo de Associação	Descrição	Interpretação Espacial
🔴 Alto-Alto (High-High)	Valor alto cercado por vizinhos também com valores altos	Indica cluster de altas magnitudes, também chamado de hotspot
🔵 Baixo-Baixo (Low-Low)	Valor baixo cercado por vizinhos com valores baixos	Indica cluster de baixas magnitudes, ou coldspot
🟠 Alto-Baixo (High-Low)	Valor alto cercado por valores baixos	Indica um outlier espacial positivo
🟡 Baixo-Alto (Low-High)	Valor baixo cercado por valores altos	Indica um outlier espacial negativo
⚪ Não significativo	Sem associação espacial relevante	O valor na região não apresenta padrão espacial detectável

Os valores do Moran Local são testados por permutação para verificar se o padrão observado é estatisticamente significativo ou poderia ocorrer por acaso.

Apenas as regiões significativas (p-valor < 0.05) costumam ser coloridas nos mapas.

Esse mapa abaixo pertence ao artigo “Are regions equal in adversity? A spatial analysis of the spread and dynamics of COVID-19 in Europe”, de Amdaoud, Arcuri e Levratto (2021). O artigo realiza uma análise espacial detalhada da mortalidade por COVID-19 em 125 regiões europeias durante a primeira onda da pandemia (março a maio de 2020).

Utilizando o índice de Moran Local (LISA), o estudo identifica padrões de autocorrelação espacial na mortalidade por COVID-19 entre as regiões europeias durante a primeira onda da pandemia. O resultado evidencia clusters espaciais significativos, com destaque para regiões High-High (altas taxas de mortalidade cercadas por outras com altas taxas), como o norte da Itália, Madrid e a região da Alsácia, na França. Ao mesmo tempo, regiões Low-Low, como o sul da Itália, Dinamarca e partes da Alemanha Oriental, apresentaram baixa mortalidade em vizinhança igualmente baixa. Esses padrões persistentes ao longo do tempo indicam que a disseminação da pandemia seguiu lógicas regionais, reforçando a importância de políticas públicas que considerem as desigualdades territoriais na resposta à crise sanitária.

Modelos de Regressão Espacial

Hipótese de independência das observações em geral é Falsa $\rightarrow$ Existe dependência espacial, ou seja, o valor observado em uma localização tende a se parecer com os valores observados em locais vizinhos.
Efeitos Espaciais $\rightarrow$ Se existir dependência ou correlação espacial significativa, devemos incluir no modelo os Efeitos Espaciais , caso contrário, as estimativas geradas pelo modelo de regressão podem ficar enviesadas, criando associações espúrias (isto é, sugerindo relações estatísticas onde elas não existem de fato, ou mascarando relações reais);
Como verificar ? $\rightarrow$ Para detectar a dependência espacial, podemos analisar os resíduos da regressão tradicional (OLS), utilizando indicadores como o índice de Moran dos resíduos. Se os resíduos mostrarem autocorrelação espacial significativa, isso indica que há estrutura espacial não modelada.
Detectada a autocorrelação espacial: e agora ? $\rightarrow$ É necessário adotar modelos de regressão que incorporem explicitamente os efeitos espaciais, evitando vieses nas estimativas dos coeficientes.

i) Modelos Globais:

Esses modelos assumem que a estrutura espacial é constante em todo o espaço geográfico. Utilizam um único parâmetro para captar a autocorrelação.
Um exemplo é o Modelo de Erro Espacial (Spatial Error Model - CAR), representado por:

\[y_i = \beta_0 + \sum^{p}_{k} \beta_k x_{ik} + \varepsilon_i \text{ , sendo } \varepsilon_i = \lambda W + \xi\]

Sendo:

$y_i$: valor observado da variável dependente na unidade $i$;
$\beta_0$: intercepto da regressão;
$\beta_k$: coeficiente da $k$-ésima variável explicativa;
$x_{ik}$: valor da $k$-ésima variável explicativa para a unidade $i$;
$\varepsilon_i$: erro com estrutura de autocorrelação espacial;
$\lambda$: coeficiente de autocorrelação espacial;
$W$: matriz de vizinhança espacial (define como cada unidade está conectada às outras);
$\xi$: erro aleatório normal (sem autocorrelação), com $

$\rightarrow$ Podemos observar que os efeitos da autocorrelação espacial são associados ao termo de erro.

🧠 Interpretação:

Se $\lambda \ne 0$, há evidência de dependência espacial nos erros, ou seja, o erro em uma localidade está correlacionado com os erros em localidades vizinhas.
Se $\lambda = 0$, o modelo reduz-se a uma regressão tradicional, e não há autocorrelação espacial significativa.

ii) Modelos Locais:

Enquanto os modelos globais assumem uma relação espacial uniforme, os modelos locais admitem que os parâmetros da regressão variam ao longo do espaço geográfico.
Um exemplo clássico é o Modelo Geograficamente Ponderado (Geographically Weighted Regression - GWR), expresso por:

\[ y_i = \beta_{0}(u_i, v_i) + \sum_{k=1}^{p} \beta_k(u_i, v_i) \cdot x_{ik} + \varepsilon_i, \\ \quad \text{com } \beta_k(u_i, v_i) \text{ estimado via } \sum_{j=1}^{n} K(d_{ij}) \cdot x_{jk} \]

Onde:

$(u_i, v_i)$: coordenadas geográficas da observação $i$;
$\beta_k(u_i, v_i)$: coeficiente local da variável explicativa $x_k$, estimado para a posição $i$;
$K(d_{ij})$: função kernel que atribui um peso à observação $j$ com base na distância $d_{ij}$ entre as localizações $i$ e $j$;

$\rightarrow$ Quanto menor $d_{ij}$, maior é o peso de $j$ na estimação de $\beta_k(u_i, v_i)$.

Fonte: ARDILLY, Pascal et al. Manuel d’analyse spatiale. 2018.

Aplicação I: Utilizando a biblioteca tmap para construção rápida de mapas temáticos

library(tmap)

# Carrega dados espaciais do mundo com variáveis socioeconômicas
data(World)

# Define o estilo do mapa (opcional, apenas visual)
tmap_style("classic")

# Cria um mapa temático utilizando a variável esperança de vida
qtm(World, fill = "life_exp", title = "Esperança de Vida por País", fill.title = "Anos",
    borders = "gray40")

Aplicação II: Baixando e construindo mapas a partir da biblioteca geobr

Bibliotecas que serão utilizadas:

library(ggplot2)
library(dplyr)
library(viridis)
library(geobr)
library(sf)
library(maptools)
library(leaflet)
library(ggspatial)

Para acessar os dados dos limites territoriais de todos os estados brasileiros é necessário utilizar a função read_state.

brasil.ufs <- read_state(code_state = "all", year = 2019, showProgress = FALSE)

Primeiro, vamos fazer um gráfico apenas com as geometrias.

ggplot(brasil.ufs) + geom_sf()

Para a construção de um mapa onde cada estado recebe uma cor de acordo com a sua região geogrática, procedemos da seguinte forma:

ggplot(brasil.ufs) + geom_sf(aes(fill = name_region)) + labs(fill = "Região")

Para plotar o mapa, agora utilizaremos dados relativos ao acesso à rede de esgoto de acordo com a unidade da federação (fonte dos dados ) segundo o censo de 2022. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.

# Entrando com os dados observados na wikipedia
acesso_esgoto <- data.frame(
  sigla = c("SP", "DF", "RJ", "MG", "ES",
            "PR", "RS", "PE", "GO", "SC",
            "BA", "MS", "RR", "CE", "SE",
            "PB", "AC", "AM", "AL", "TO",
            "MT", "RN", "PI", "PA", "RO",
            "AP", "MA"),
  
  rede_esgoto = c(93.6, 89.9, 88.9, 81.4, 78.2,
                  73.7, 71.0, 65.1, 61.7, 61.4,
                  59.5, 57.2, 56.0, 52.8, 51.0,
                  49.9, 47.7, 47.0, 40.6, 39.8,
                  37.4, 35.8, 26.8, 23.5, 21.7,
                  20.7, 18.1)
)

# Unindo os dados com o shapefile

brasil.ufs <- brasil.ufs %>%
    left_join(acesso_esgoto, by = c(abbrev_state = "sigla"))

# Construindo o mapa com ggplot
brasil.ufs |>
    ggplot(aes(fill = rede_esgoto), color = "black") + geom_sf() + scale_fill_viridis(name = "Municípios com rede de esgoto (%)",
    direction = -1) + xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") +
    annotation_north_arrow(location = "br") + theme_minimal()

Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF).

# Calcula o centroide de cada estado no objeto espacial brasil.ufs.
coord_pontos <- brasil.ufs |>
    st_centroid()

# Construindo o mapa com ggplot
ggplot(brasil.ufs) + geom_sf() + geom_sf(data = coord_pontos, aes(size = rede_esgoto),
    col = "blue", alpha = 0.65, show.legend = "point") + scale_size_continuous(name = "Municípios com rede de esgoto (%)") +
    xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") + annotation_north_arrow(location = "br") +
    theme_minimal()

Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet

leaflet(coord_pontos) |>
  # Camada tradicional com nomes e divisões políticas
  addProviderTiles(providers$OpenStreetMap, group = "Mapa Político") |>
  # Camada de imagem de satélite
  addProviderTiles(providers$Esri.WorldImagery, group = "Satélite") |>
  # Adiciona os círculos nos centroides
  addCircleMarkers(
    label = ~ paste0(abbrev_state, ": ", rede_esgoto, "%"), # Rótulo exibido ao passar o mouse
    labelOptions = labelOptions(textsize = "13px"), # Tamanho da fonte do rótulo
    radius = ~ rede_esgoto / 10,              # Raio do círculo proporcional
    fillOpacity = 0.5,                       # Transparência dos círculos
    group = "Dados"
  ) |>
  # Adiciona controle para alternar entre as camadas de fundo
  addLayersControl(
    baseGroups = c("Mapa Político", "Satélite"),
    overlayGroups = "Dados",
    options = layersControlOptions(collapsed = FALSE)
  )

Aplicação III: Dengue em Dourados/MS - Parte 1: Análise exploratória

Nesta apresentação, serão utilizados dados com ruído espacial, com o objetivo de preservar a confidencialidade das informações originais. Esses dados fazem parte da dissertação de Isis Rodrigues Reitman, intitulada “Saúde e Ambiente Urbano: a relação de incidência de dengue e as disparidades espaciais em Dourados – MS”, apresentada em abril de 2016 no Programa de Pós-Graduação em Geografia da Universidade Federal da Grande Dourados (UFGD).
Os ruídos nas coordenadas foram gerados a partir de uma distribuição uniforme. A função produz valores aleatórios entre -15 e 15 para cada linha do conjunto de dados. Esse procedimento simula um deslocamento máximo de até 15 metros em cada direção, o que é suficiente para preservar o padrão espacial geral dos casos, ao mesmo tempo em que evita a identificação precisa das localizações originais.

# Supondo que seu dataframe se chame `casos` Adicionando um ruído aleatório
# pequeno (por exemplo, até ±15 metros) às coordenadas X e Y

set.seed(123)  # Para reprodutibilidade
casos$X_ruido <- casos$X + runif(nrow(casos), -15, 15)
casos$Y_ruido <- casos$Y + runif(nrow(casos), -15, 15)

Essa aplicação também se encontra no Curso de Estudos Ecológicos ministrado para o curso de Pós-Graduação em Epidemiologia em Saúde Pública em 2019, pelos pesquisadores Oswaldo Gonçalves Cruz (PROCC/FIOCRUZ) e Wagner Tassinari (DEMAT/ICE/UFRRJ)

OBS: Os dados e as malhas geográficas utilizadas nessa apresentação, estão disponíveis no seguinte endereço: (link)

Biliotecas do R que serão utilizadas

library(readr)
library(tidyverse)
library(sf)
library(maptools)
library(spatstat)
library(tmap)

Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS

pop2010 <- read_csv("dados/dengue_dourados/pop2010.csv")
setor.sf <- read_sf("dados/dengue_dourados/Setor_UTM_SIRGAS.shp", crs = 31981)
contorno.sf <- read_sf("dados/dengue_dourados/contorno.shp", crs = 31981)

📌 O que é CRS (Coordinate Reference System) ?

O CRS é como uma “regra de tradução” entre o que está em um mapa e o mundo real. Ele define como as coordenadas (como 140, 12) se relacionam com locais verdadeiros na Terra, dizendo se os valores estão em metros, graus ou outra unidade, e onde fica o ponto de partida (a origem). Sem um CRS, uma coordenada não tem significado, pois não sabemos onde ela realmente está nem em que escala.

🗺️ Por que isso é importante ?

Se dois mapas tiverem CRSs diferentes, os pontos não vão se alinhar corretamente — como tentar juntar peças de quebra-cabeça de caixas diferentes. Usar o CRS certo garante que todas as informações espaciais estejam bem posicionadas e coerentes. Por exemplo, o WGS 84 é usado por GPS e Google Maps (em graus), enquanto o UTM é usado para mapas locais mais precisos (em metros). Saber isso evita erros em análises espaciais e garante que os dados “conversem” entre si.

OBS: Se quiser saber um pouco mais a respeito do CRS, basta acessar (link)

Fazendo um join com as tabelas com os setores censitários + população

setor.sf <- setor.sf |>
    mutate(idsetor = as.numeric(CD_GEOCODI)) |>
    inner_join(pop2010, by = "idsetor")

Lendo e plotando os casos de dengue georreferenciados em Dourados/MS

casos <- read_csv("dados/dengue_dourados/dengue_dourados.csv")
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = casos.sf,
    color = "red", size = 1) + theme_void()

Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS

lixo <- read_csv2("dados/dengue_dourados/lixo_dourados.csv")
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo.sf,
    color = "blue", size = 1) + theme_void()

Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS

Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando `st_intersection`.

lixo2.sf <- st_intersection(contorno.sf, lixo.sf)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo2.sf,
    color = "blue", size = 1) + theme_void()

Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)

# Adiciona categorias para legenda
lixo$tipo <- "Lixo"
casos.sf$tipo <- "Caso"

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = casos.sf, color = "red",
    size = 0.7, aes(colour = "Caso"), show.legend = "point") + geom_sf(data = lixo2.sf,
    color = "salmon", size = 1, aes(colour = "Lixo"), show.legend = "point") + scale_fill_distiller(palette = "PuBu",
    direction = 1) + scale_colour_manual(values = c(Caso = "red", Lixo = "salmon")) +
    theme_minimal()

### Construindo os buffers

Serão construídos buffers com raio de 500 metros ao redor de cada ponto de coleta de lixo. Essa abordagem permite analisar se os casos de dengue ocorrem dentro desse perímetro, ou seja, a até 500 metros dos pontos de coleta. O objetivo é investigar se existe uma relação entre a proximidade desses locais e a ocorrência de casos de dengue.

Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.

lixo_buffer <- st_buffer(lixo2.sf, 500)

# Adiciona categorias para legenda
lixo_buffer$tipo <- "Lixo"
casos.sf$tipo <- "Caso"

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = lixo_buffer, aes(color = tipo),
    fill = "transparent", size = 0.4) + geom_sf(data = casos.sf, aes(color = tipo),
    size = 0.7) + scale_fill_distiller(palette = "PuBu", direction = 1) + scale_color_manual(values = c(Caso = "red",
    Lixo = "gray")) + theme_minimal()

Representando os casos e o lixo de forma interativa com o a biblioteca `tmap`

library(tmap)


tmap_mode("view")  # modo interativo (Leaflet)
# tmap_mode('plot') # modo estático (p/ publicação)

# Constrói o mapa
tm_shape(setor.sf) + tm_borders("black") + tm_shape(casos.sf) + tm_dots(col = "red",
    palette = "red") + tm_shape(lixo.sf) + tm_dots(col = "green", palette = "green") +
    tm_shape(lixo_buffer) + tm_borders(col = "blue")

Convertendo o dado de pontos (padrão pontual) para dados de área

library(sf)
library(dplyr)

# Conta casos confirmados (CLASSI_FIN == 1) por setor
casos <- casos.sf |>
  filter(CLASSI_FIN == 1) |>        # Filtra apenas casos confirmados
  mutate(contador = 1) |>          # Adiciona coluna de contagem
  st_join(setor.sf, join = st_within) |>  # Une os casos aos setores (dentro do polígono)
  group_by(ID) |>                 # Agrupa por ID do setor
  summarise(casos = sum(contador, na.rm = TRUE)) |> # Conta casos por setor
  st_drop_geometry()  ## remove atributos de geometria


# Número de depósitos de lixo por setor (sem geometria)
lixo <- lixo.sf |>
  mutate(contador = 1) |>                      # Adiciona coluna de contagem
  st_join(setor.sf, join = st_within) |>       # Junta pontos aos setores onde estão contidos
  group_by(ID) |>                              # Agrupa por setor
  summarise(lixo = sum(contador, na.rm = TRUE)) |>  # Soma os depósitos por setor
  st_drop_geometry()                           # Remove geometria


# Inserindo as contagens de lixo e casos confirmados na geometria dos setores
setor.sf <- setor.sf |>
  left_join(lixo, by = "ID") |>    # Junta a contagem de depósitos de lixo por setor
  left_join(casos, by = "ID")      # Junta a contagem de casos confirmados por setor

Plotando o mapa temático dos casos por setor censitário

plot(setor.sf["casos"])

Plotando o mapa temático dos pontos de coleta de lixo por setor censitário

plot(setor.sf["lixo"])

Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário

setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0  # Transformando os missings em zero

summary(setor.sf$tx)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.784   4.075   4.527  56.399

Plotando a distribuição da incidência em Dourados/MS

library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")

ggplot(setor.sf) + geom_sf(aes(fill = tx), color = "black") + scale_fill_gradientn(colours = pal) +
    ggtitle("Taxa de incidência de Dengue") + theme_void()

Kernel por atributos

Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.
Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA

dourados.contorno <- st_union(setor.sf)
plot(dourados.contorno)

- Agora iremos converter o contorno geográfico da cidade de Dourados (um polígono sf) para o formato owinda classe spatstat, que é usado para análise espacial de padrões de pontos.

dourados.w <- as.owin(st_geometry(dourados.contorno))

Extraindo os centróides dos polígonos em Dourados/MS

centroides <- st_centroid(st_geometry(setor.sf))

# Transformando em os centróides em formato sp (spatstat)
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c("X", "Y")

plot(centroides.sp)

Colocando os pontos no formato sp

centroides.ppp <- ppp(centroides.sp$X, centroides.sp$Y, dourados.w)

plot(centroides.ppp, pch = 19, cex = 0.5)

Construindo o kernel por atributo da taxa de detecção

kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)

centroides.ppp: Objeto do tipo ppp com os pontos centrais (centroides) dos setores. 500 Largura de banda (bandwidth), ou seja, o alcance de suavização da função kernel, neste caso 500 metros (ou unidades do sistema).
weights = setor.sf$tx: Peso de cada ponto (setor), neste caso uma variável chamada tx, que provavelmente representa uma taxa (ex: incidência por população).
scalekernel = TRUE: Faz a função kernel se ajustar automaticamente à soma total dos pesos, preservando a escala da densidade ponderada.

Construindo a matriz de vizinhança para verificar a autocorrelação espacial

library(spdep)
viz <- poly2nb(setor.sf)
viz

Neighbour list object:
Number of regions: 284 
Number of nonzero links: 1728 
Percentage nonzero weights: 2.142432 
Average number of links: 6.084507

setor.sf é um objeto do tipo sf com polígonos (setores).
poly2nb() (polygon to neighbors) calcula quais polígonos são vizinhos, ou seja, compartilham ao menos um lado ou ponto de fronteira (por padrão).
viz é um objeto do tipo nb (“neighbors list”) — uma lista onde cada item contém os índices dos setores vizinhos ao setor correspondente.

Iremos precisar da coordenadas dos centróides para montar a malha de conectividade.

setor.sp <- as(setor.sf, "Spatial")  # convertendo em formato sp
coord <- coordinates(setor.sp)  # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

Verificando a malha de conectividade da vizinhança de Dourados/MS

viz.sf <- as(nb2lines(viz, coords = coord), "sf")
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))

st_set_crs() define o sistema de referência de coordenadas (CRS) do objeto viz.sf.
Ele está copiando o CRS do objeto setor.sf, garantindo que os dois objetos estejam no mesmo sistema de coordenadas.

# Plotando o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) + geom_sf(fill = "lightpink", color = "white") + geom_sf(data = viz.sf) +
    theme_minimal() + ggtitle("Vizinhança por \n conectividade") + ylab("Latitude") +
    xlab("Longitude")
mapa.viz

Obtendo a correlação da taxa de incidência de dengue Dourados/MS

Convertendo a estrutura criada por poly2nb(setor.sf) um formato que funções como moran.test.

pesos.viz <- nb2listw(viz)
moran.test(setor.sf$tx, pesos.viz)


    Moran I test under randomisation

data:  setor.sf$tx  
weights: pesos.viz    

Moran I statistic standard deviate = 16.09, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.536654789      -0.003533569       0.001127179

Podemos observar que existe uma forte autocorrelação espacial positiva na variável tx, ou seja, valores semelhantes estão espacialmente agrupados. Em outras palavras, setores com valores altos de tx tendem a estar próximos de outros setores com valores igualmente altos (e o mesmo vale para valores baixos). Esse padrão não é aleatório, já que o p-valor é muito pequeno, rejeitando a hipótese nula de ausência de autocorrelação espacial.

Plotando o correlograma

O correlograma espacial de Moran’s I que você gerou está te mostrando como a autocorrelação espacial da variável tx varia conforme o “distanciamento” espacial entre os setores.

correl <- sp.correlogram(viz, setor.sf$tx, order = 8, method = "I")
correl

Spatial correlogram for setor.sf$tx 
method: Moran's I
           estimate expectation    variance standard deviate Pr(I) two sided
1 (284)  0.53665479 -0.00353357  0.00112718          16.0897       < 2.2e-16
2 (284)  0.25052294 -0.00353357  0.00048462          11.5407       < 2.2e-16
3 (284)  0.10453394 -0.00353357  0.00028323           6.4213       1.351e-10
4 (284) -0.01860594 -0.00353357  0.00019717          -1.0734         0.28309
5 (284) -0.07828447 -0.00353357  0.00016719          -5.7811       7.423e-09
6 (284) -0.06556547 -0.00353357  0.00016167          -4.8787       1.068e-06
7 (284) -0.03219851 -0.00353357  0.00016741          -2.2154         0.02673
8 (284) -0.02957621 -0.00353357  0.00018342          -1.9229         0.05449
           
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)    
5 (284) ***
6 (284) ***
7 (284) *  
8 (284) .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

-A tabela consiste:

estimate: o valor observado de Moran’s I.
expectation: valor esperado sob aleatoriedade (≈ 0).
Pr(I): p-valor de significância para cada lag.

plot(correl)

Nesse gráfico podemos observar:
1. O eixo y mostra o valor de Moran’s I para cada lag.
2. As barras de erro indicam o intervalo de confiança do Moran’s I em cada lag.
Podemos compluir então que:
- Lag 1 (vizinhos imediatos): Moran’s I = 0.5367, altamente significativo (p < 0.001) $\rightarrow$ forte autocorrelação espacial positiva.
- Lags 2 e 3: valores positivos e ainda significativos $\rightarrow$ padrão espacial persiste até o terceiro nível de vizinhança.
- Lags 4 a 8: valores negativos ou próximos de zero, com significância decrescente ou ausência de significância $\rightarrow$ a autocorrelação desaparece ou se inverte à medida que a distância aumenta.

Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.

setor.sf$pval <- localmoran(setor.sf$tx, pesos.viz)[, 5]

tm_shape(setor.sf) + tm_polygons(col = "pval", title = "p-valores", breaks = c(0,
    0.01, 0.05, 0.1, 1), style = "fixed", palette = "-Oranges", border.col = "grey")

Moran Local (Lisa Map) da taxa de incidência Dourados/MS

resI <- localmoran.sad(lm(setor.sf$tx ~ 1), 1:length(viz), viz, style = "W")
summary(resI)[1:10, ]

    Local Morans I Stand. dev. (N)   Pr. (N) Saddlepoint Pr. (Sad)  Expectation
1 1    -0.01015882     -0.01643092 1.0131094 -0.03718329 0.9703389 -0.003533569
2 2     0.19345824      0.39747581 0.6910166  0.61571046 0.5380856 -0.003533569
3 3    -0.01344306     -0.02457595 1.0196068 -0.05129549 0.9590901 -0.003533569
4 4     0.53333749      1.44072903 0.1496612  1.52458273 0.1273632 -0.003533569
5 5     0.00291006      0.01729191 0.9862037  0.02065198 0.9835233 -0.003533569
     Variance    Skewness Kurtosis   Minimum  Maximum         omega       sad.r
1 1 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001428817 -0.01636976
2 2 0.2456263 -0.04188294 8.752890 -70.87400 69.87400  0.0024931754  0.38392781
3 3 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0002135955 -0.02448076
4 4 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0068568003  1.16560992
5 5 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0001628523  0.01723265
          sad.u
1 1 -0.01637534
2 2  0.41965894
3 3 -0.02449684
4 4  1.77121677
5 5  0.01723367
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]

A função localmoran.sad() calcula o Moran Local (Iᵢ) e algumas outras estatísticas associadas para cada polígono (setor).
lm(setor.sf$tx ~ 1): cria um modelo linear sem preditores (apenas o intercepto), ou seja, considera a variável tx como resposta a ser analisada.
1:length(viz): identifica o índice de cada área/setor.
viz: é a estrutura de vizinhança (do tipo nb), que define quem é vizinho de quem.
style = "W": define o estilo de normalização dos pesos espaciais (padronização row-standardized).

A está mostrando:
- Local Morans I: Valor do I de Moran local (Iᵢ), indicando se o valor do setor e seus vizinhos são semelhantes (positivo) ou diferentes (negativo).
- Stand. dev. (N): Desvio padrão padronizado da estatística de Moran local.
- Pr. (N): p-valor do teste de significância para o valor de Iᵢ (baseado na dist. normal).
- Saddlepoint: Valor do I de Moran local corrigido.

setor.sf$MoranLocal <- summary(resI)[, 1]

library(scales)

ggplot(setor.sf) + geom_sf(aes(fill = MoranLocal), color = "black") + scale_fill_gradientn(colours = c("blue",
    "white", "red"), values = rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))),
    guide = "colorbar") + ggtitle("Moran local") + theme_void()

Neste mapa, podemos observar:

Valores altos de Local Morans I com p-valor pequeno indicam clusters estatisticamente significativos:
- Positivo + significativo $\rightarrow$ cluster alto-alto ou baixo-baixo.
- Negativo + significativo $\rightarrow$ outlier (ex: um valor alto cercado por baixos, ou o contrário).

LISA Map

library(spdep)

# Convertendo uma estrutura de vizinhança (objeto nb) em uma matriz de pesos
# espaciais padronizada
w <- nb2listw(viz, style = "W")

# Matriz de pesos espaciais padronizada (row-standardized)
resI <- localmoran(setor.sf$tx, nb2listw(viz, style = "W"))

# Armazena estatísticas no objeto espacial
setor.sf$Ii <- resI[, 1]  # valor de Moran Local
setor.sf$pvalue <- resI[, 5]  # p-valor

# Valor padronizado da variável tx
setor.sf$z_tx <- scale(setor.sf$tx)[, 1]

# Média dos vizinhos (lagged value)
lag_tx <- lag.listw(w, setor.sf$z_tx)

# Classificação dos tipos de associação local dos clusters
setor.sf$cluster_type <- "Não significativo"

setor.sf$cluster_type[setor.sf$pvalue <= 0.05 & setor.sf$z_tx > 0 & lag_tx > 0] <- "Alto-Alto"
setor.sf$cluster_type[setor.sf$pvalue <= 0.05 & setor.sf$z_tx < 0 & lag_tx < 0] <- "Baixo-Baixo"
setor.sf$cluster_type[setor.sf$pvalue <= 0.05 & setor.sf$z_tx > 0 & lag_tx < 0] <- "Alto-Baixo"
setor.sf$cluster_type[setor.sf$pvalue <= 0.05 & setor.sf$z_tx < 0 & lag_tx > 0] <- "Baixo-Alto"

library(tmap)

# Modo estático (ideal para relatórios impressos ou PDF)
tmap_mode("plot")

# Tema visual limpo com fundo branco
tmap_style("white")

# Geração do mapa LISA
tm_shape(setor.sf) +
  tm_fill(
    col = "cluster_type",
    palette = c(
      "Alto-Alto"         = "#E60000",  # vermelho forte
      "Baixo-Baixo"       = "#0033CC",  # azul escuro
      "Baixo-Alto"        = "#9999FF",  # azul claro
      "Alto-Baixo"        = "#FF9999",  # rosa claro
      "Não significativo" = "#FFFFFF"   # branco
    ),
    title = "Cluster Local (LISA)",
    legend.is.portrait = TRUE
  ) +
  tm_borders(col = "gray40", lwd = 0.4) +
  tm_layout(
    main.title = "Mapa LISA - Clusters Locais da Taxa de Incidência",
    main.title.size = 1.2,
    legend.outside = TRUE,
    frame = FALSE,
    bg.color = "white"
  )

Aplicação IV: Dengue em Dourados/MS - Parte 2: Modelagem (Modelos Linear, CAR e GWR)

Ajustando o modelo de regressão linear simples.

# Transformando os missings em zero
setor.sf$lixo[is.na(setor.sf$lixo)] <- 0

dourados.lm <- lm(tx ~ lixo, data = setor.sf)
summary(dourados.lm)


Call:
lm(formula = tx ~ lixo, data = setor.sf)

Residuals:
   Min     1Q Median     3Q    Max 
-4.492 -4.122 -2.272  0.262 51.907 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.4917     0.4779   9.400   <2e-16 ***
lixo         -0.3697     0.1901  -1.944   0.0528 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.198 on 282 degrees of freedom
Multiple R-squared:  0.01323,   Adjusted R-squared:  0.009731 
F-statistic: 3.781 on 1 and 282 DF,  p-value: 0.05283

Checando os residuos para verificar a presença de autocorrelação.

dourados.lm$lmresid <- residuals(dourados.lm)
moran.test(dourados.lm$lmresid, pesos.viz)


    Moran I test under randomisation

data:  dourados.lm$lmresid  
weights: pesos.viz    

Moran I statistic standard deviate = 15.949, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.532194688      -0.003533569       0.001128294

Ajustando o modelo CAR (Spatial Error Model)

library(spatialreg)
dourados.car <- errorsarlm(tx ~ lixo, data = setor.sf, listw = pesos.viz)
summary(dourados.car)


Call:errorsarlm(formula = tx ~ lixo, data = setor.sf, listw = pesos.viz)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.2598  -2.4487  -1.0865   1.1696  32.1944 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.99554    1.56936  2.5460   0.0109
lixo        -0.21893    0.14217 -1.5399   0.1236

Lambda: 0.81645, LR test value: 176.59, p-value: < 2.22e-16
Asymptotic standard error: 0.040221
    z-value: 20.299, p-value: < 2.22e-16
Wald statistic: 412.05, p-value: < 2.22e-16

Log likelihood: -874.2534 for error model
ML residual variance (sigma squared): 23.437, (sigma: 4.8412)
Number of observations: 284 
Number of parameters estimated: 4 
AIC: 1756.5, (AIC for lm: 1931.1)

Checando os residuos para verificar a presença de autocorrelação

dourados.car$carresid <- residuals(dourados.car)
moran.test(dourados.car$carresid, pesos.viz)


    Moran I test under randomisation

data:  dourados.car$carresid  
weights: pesos.viz    

Moran I statistic standard deviate = 0.66469, p-value = 0.2531
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.018969292      -0.003533569       0.001146147

Ajustando o modelo GWR (Geographically Weighted Regression)

# Biblioteca para ajustar o modelos GWR
library(spgwr)

# Estimando a largura de banda “ideal” para o kernel
GWRbanda <- gwr.sel(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y),
    adapt = T)

Adaptive q: 0.381966 CV score: 14008.03 
Adaptive q: 0.618034 CV score: 14380.73 
Adaptive q: 0.236068 CV score: 13508.92 
Adaptive q: 0.145898 CV score: 12649.04 
Adaptive q: 0.09016994 CV score: 11470.56 
Adaptive q: 0.05572809 CV score: 10305.47 
Adaptive q: 0.03444185 CV score: 9396.074 
Adaptive q: 0.02128624 CV score: 9035.635 
Adaptive q: 0.01315562 CV score: 8561.705 
Adaptive q: 0.008130619 CV score: 8237.524 
Adaptive q: 0.005024999 CV score: 8159.145 
Adaptive q: 0.003965528 CV score: 8248.869 
Adaptive q: 0.006099679 CV score: 8196.208 
Adaptive q: 0.00543549 CV score: 8163.423 
Adaptive q: 0.005086408 CV score: 8158.764 
Adaptive q: 0.005127099 CV score: 8158.728 
Adaptive q: 0.005167789 CV score: 8158.857 
Adaptive q: 0.005127099 CV score: 8158.728

# Ajustando o modelo GWR
dourados.gwr = gwr(tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, centroides.sp$Y),
    adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)

dourados.gwr

Call:
gwr(formula = tx ~ lixo, data = setor.sf, coords = cbind(centroides.sp$X, 
    centroides.sp$Y), adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)
Kernel function: gwr.Gauss 
Adaptive quantile: 0.005127099 (about 1 of 284 data points)
Summary of GWR coefficient estimates at data points:
                    Min.     1st Qu.      Median     3rd Qu.        Max.
X.Intercept.  -0.0656658   1.4352539   2.9363259   4.8308391  38.4197448
lixo         -20.4306144  -1.2551754  -0.3598796  -0.0063376   1.5698378
              Global
X.Intercept.  4.4917
lixo         -0.3697
Number of data points: 284 
Effective number of parameters (residual: 2traceS - traceS'S): 139.5371 
Effective degrees of freedom (residual: 2traceS - traceS'S): 144.4629 
Sigma (residual: 2traceS - traceS'S): 5.129984 
Effective number of parameters (model: traceS): 104.7659 
Effective degrees of freedom (model: traceS): 179.2341 
Sigma (model: traceS): 4.605575 
Sigma (ML): 3.658769 
AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 1881.684 
AIC (GWR p. 96, eq. 4.22): 1647.491 
Residual sum of squares: 3801.792 
Quasi-global R2: 0.7432628

# Colocando a saída do modelo dentro de um objeto dataframe.
results <- as.data.frame(dourados.gwr$SDF)
head(results)

     sum.w X.Intercept.       lixo X.Intercept._se  lixo_se       gwr.e
1 4.362973     6.333142 -1.7875262        2.310328 1.706014 -0.08071313
2 3.861387     7.610174 -0.9533830        2.360779 1.550726  2.25824677
3 3.059017     6.825286 -1.0059712        2.605365 1.524275 -1.97316059
4 3.867575     9.254170 -0.4391772        2.287048 1.215040  5.07389644
5 5.343561     4.415895 -0.5911405        2.008245 1.129649 -2.55248940
      pred  pred.se   localR2 X.Intercept._se_EDF lixo_se_EDF pred.se.1
1 2.758089 2.466160 0.6473012            2.573392    1.900268  2.746967
2 7.610174 2.360779 0.5332692            2.629587    1.727299  2.629587
3 5.819314 2.028659 0.5115313            2.902022    1.697835  2.259651
4 8.814992 1.742409 0.3882650            2.547460    1.353390  1.940807
5 3.824754 1.458803 0.6828771            2.236912    1.258276  1.624909
   coord.x coord.y
1 731406.5 7541547
2 730845.4 7541481
3 730957.7 7541239
4 730487.2 7541437
5 729874.1 7541446
 [ reached 'max' / getOption("max.print") -- omitted 1 rows ]

Verificando a distribuição dos coeficientes de regressão para a variável lixo

hist(results$lixo)
abline(v = median(results$lixo), col = "red")

Verificando a distribuição dos localR2

hist(results$localR2)
abline(v = median(results$localR2), col = "blue")

Incorporando alguns parâmetros de saída do modelo na tabela setor.sf

setor.sf$coef.lixo <- results$lixo
setor.sf$localR2 <- results$localR2
setor.sf$pred.gwr <- results$pred

Mapa dos coeficientes de regressão para a variável lixo

map.lixo <- ggplot(setor.sf) + geom_sf(aes(fill = coef.lixo), color = "black") +
    scale_fill_gradientn(colours = pal) + ggtitle("Distribuição dos coef. var. lixo") +
    theme_void()
map.lixo

Checando os residuos para verificar a presença de autocorrelação para o modelo GWR.

# Calculando os resíduos para o modelo GWR
results$residuos <- setor.sf$tx - results$pred

moran.test(results$residuos, pesos.viz)


    Moran I test under randomisation

data:  results$residuos  
weights: pesos.viz    

Moran I statistic standard deviate = 1.4555, p-value = 0.07277
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.046011652      -0.003533569       0.001158754

Mapeando os coeficientes de regressão para a variável lixo por significancia através do teste de wald

# Calculando a estatística de wald
setor.sf$wald.teste <- abs(results$lixo/results$lixo_se)
# Dicotomizando a estatística de wald
setor.sf$wald.teste <- ifelse(setor.sf$wald.teste < 2, 0, 1)


map.wald <- ggplot(setor.sf) + geom_sf(aes(fill = factor(wald.teste)), color = "black") +
    scale_fill_manual(values = c("white", "purple"), labels = c("< 2", ">= 2"), name = "Wald") +
    ggtitle("Coef. lixo significativos") + theme_void()


library(gridExtra)
grid.arrange(map.lixo, map.wald, ncol = 2)

Mapa dos coeficientes de determinação regionalizados ($R^2$ local).

ggplot(setor.sf) + geom_sf(aes(fill = localR2), color = "black") + scale_fill_gradientn(colours = pal) +
    ggtitle("R² local") + theme_void()

Verificando a distribuição dos preditos.

histdens <- function(x, titulo = "") {
    densi <- density(x)
    xli <- range(densi$x)
    yli <- range(densi$y)
    hist(x, col = "red", probability = T, xlim = xli, ylim = yli, main = titulo)
    lines(densi, lwd = 2)
    abline(v = median(x), lwd = 2, col = 4, lty = 2)
}

par(mfrow = c(2, 2))

hist.tx <- histdens(setor.sf$tx, "Tx Bruta")
hist.lm <- histdens(dourados.lm$fitted.values, "Pred LM")
hist.car <- histdens(dourados.car$fitted.values, "Pred CAR")
hist.gwr <- histdens(results$pred, "Pred GWR")

Mapeando os valores observados e preditos dos modelos ajustados

library(colorspace)  # 

setor.sf$brks <- cut(setor.sf$tx, include.lowest = TRUE, right = TRUE, breaks = c(-4,
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))

tx.bruta <- ggplot(setor.sf) + geom_sf(aes(fill = brks), color = "black") + ggtitle("Taxa Bruta") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()



setor.sf$brks.lm <- cut(dourados.lm$fitted.values, lowest = TRUE, right = TRUE, breaks = c(-4,
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.lm), color = "black") + ggtitle("Taxa Predita - LM") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()

setor.sf$brks.car <- cut(dourados.car$fitted.values, lowest = TRUE, right = TRUE,
    breaks = c(-4, 0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10",
        "> 10"))


pred.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - CAR") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()

setor.sf$brks.gwr <- cut(results$pred, lowest = TRUE, right = TRUE, breaks = c(-4,
    0, 2, 4, 10, 57), labels = c("0", "0 - 2", "2 - 4", "4 - 10", "> 10"))


pred.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.car), color = "black") + ggtitle("Taxa Predita - GWR") +
    scale_fill_discrete_sequential(palette = "Heat2", c1 = 80, c2 = 30, l1 = 30,
        l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75", drop = FALSE, name = "Taxa") +
    theme_void()

library(gridExtra)
grid.arrange(tx.bruta, pred.lm, pred.car, pred.gwr, ncol = 2)

Comparando a distribuição dos resíduos dos modelos ajustados

library(vioplot)
vioplot(dourados.lm$residuals, dourados.car$residuals, results$residuos, names = c("LM",
    "CAR", "GWR"), col = "orange")
title("Gráficos violinos da distribuição dos resíduos")

Mapeando os resíduos dos modelos ajustados

library(colorspace)  # 

setor.sf$brks.res.lm <- cut(dourados.lm$residuals, include.lowest = TRUE, right = TRUE,
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
        "> 5"))

res.lm <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.lm), color = "black") +
    ggtitle("Resíduos - LM") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
    drop = FALSE, name = "Taxa") + theme_void()

setor.sf$brks.res.car <- cut(dourados.car$residuals, include.lowest = TRUE, right = TRUE,
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
        "> 5"))

res.car <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.car), color = "black") +
    ggtitle("Resíduos - CAR") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
    drop = FALSE, name = "Taxa") + theme_void()


setor.sf$brks.res.gwr <- cut(results$residuos, include.lowest = TRUE, right = TRUE,
    breaks = c(-14, -5, -1, 1, 5, 52), labels = c("< -5", "-5 a -1", "0", "1 a 5",
        "> 5"))

res.gwr <- ggplot(setor.sf) + geom_sf(aes(fill = brks.res.gwr), color = "black") +
    ggtitle("Resíduos - GWR") + scale_fill_discrete_sequential(palette = "Purple-Yellow",
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5, na.value = "grey75",
    drop = FALSE, name = "Taxa") + theme_void()


library(gridExtra)
grid.arrange(res.lm, res.car, res.gwr, ncol = 2)

Aplicação V: Geoestatística com dados de pluviosidade na cidade do Rio de Janeiro/RJ

Biliotecas do R que serão utilizadas

library(sf)
library(sp)
library(dplyr)
library(gstat)
library(lattice)
library(automap)
library(raster)
library(leaflet)

Importando a tabela com a chuva acumulada média de 7 dias das últimas 24hs e 96hs das estações pluviométricas da cidde do Rio de Janeiro

Descrição das Estações (Alerta Rio)

Download Dados Pluviométricos (Alerta Rio)

pluvio <- read.csv("dados/chuva_rio/pluviosidade.csv", sep = ";")

Análise Gráfica Descritiva

quantil <- quantile(pluvio$acumulado_24h, seq(0, 1, 0.2))
quantil

   0%   20%   40%   60%   80%  100% 
 0.00  1.92  4.24  7.36 11.84 36.00

Transformar os dados em um objeto espacial do R

x - Longitude
y - Latitude

sp::coordinates(pluvio) = ~x + y

Análise Gráfica Descritiva com os dados espaciais

# Bubble plot
bubble(pluvio, "acumulado_24h", key.entries = quantil, pch = 19, col = "blue")

# Point plot
spplot(pluvio["acumulado_24h"], scales = list(draw = T), key.space = "right", colorkey = T)

Modelando o variograma experimental (ou empírico)

width - Distância média entre amostras ou distância dos lags
cutoff - Máxima distância

variogram.emp = variogram(acumulado_24h ~ x + y, pluvio, width = 1000, cutoff = 20000)
variogram.emp

   np      dist     gamma dir.hor dir.ver   id
1   2  1385.875 20.381155       0       0 var1
2   6  2543.894  8.139327       0       0 var1
3   8  3648.158  5.261490       0       0 var1
4  15  4583.562 13.707344       0       0 var1
5  22  5491.324 26.228076       0       0 var1
6  15  6531.432 52.574023       0       0 var1
7  12  7512.398 42.179960       0       0 var1
8  20  8516.447 65.283855       0       0 var1
9  15  9432.920 36.022491       0       0 var1
10 19 10450.957 36.569808       0       0 var1
11 13 11492.778 76.150012       0       0 var1
12 19 12568.517 39.354425       0       0 var1
 [ reached 'max' / getOption("max.print") -- omitted 7 rows ]

# Variogram plot
plot(variogram.emp, main = "Empirical variogram", pch = 19, col = "darkblue")

Ajustando o semivariograma teórico

Sill - Semivariância estrutural ou contribuição (ponto máximo que chega ao plato no eixo de y)
Range - Alcance (ponto máximo em x)
Nugget - Efeito pepita, quando (y, 0)

Fonte do gráfico

## 1) Modelo Linear
lin.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Lin", range = 5000,
    nugget = 10))

## 2) Modelo exponencial
exp.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Exp", range = 5000,
    nugget = 10))

## 3) Modelo gaussiano
gau.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Gau", range = 5000,
    nugget = 10))

## 3) Modelo wave
wav.fit = fit.variogram(variogram.emp, model = vgm(psill = 240, model = "Wav", range = 5000,
    nugget = 10))

plot.lin <- plot(variogram.emp, lin.fit, main = "Modelo Linear")
plot.exp <- plot(variogram.emp, exp.fit, main = "Modelo Exponencial")
plot.gau <- plot(variogram.emp, gau.fit, main = "Modelo Gaussiano")
plot.wav <- plot(variogram.emp, wav.fit, main = "Modelo Wave")

Validação cruzada

RMSE (root mean squared error): é a medida que calcula “a raiz quadrática média” dos erros entre valores observados (reais) e predições (hipóteses).

## 1) Modelo Linear
cv.lin <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = lin.fit)
summary(cv.lin)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred          var1.var         observed        residual        
 Min.   : 0.6757   Min.   : 9.001   Min.   : 0.00   Min.   :-10.41265  
 1st Qu.: 4.3063   1st Qu.:11.956   1st Qu.: 2.20   1st Qu.: -2.84581  
 Median : 6.7569   Median :19.120   Median : 5.00   Median : -0.44588  
 Mean   : 7.2750   Mean   :28.832   Mean   : 7.34   Mean   :  0.06504  
 3rd Qu.:10.4716   3rd Qu.:43.382   3rd Qu.:10.30   3rd Qu.:  1.60176  
 Max.   :15.6127   Max.   :80.832   Max.   :36.00   Max.   : 29.36110  
     zscore                fold      
 Min.   :-2.0312467   Min.   : 1.00  
 1st Qu.:-0.5694160   1st Qu.: 8.25  
 Median :-0.1237453   Median :15.50  
 Mean   : 0.0004815   Mean   :15.50  
 3rd Qu.: 0.3209767   3rd Qu.:22.75  
 Max.   : 4.1304203   Max.   :30.00

plot(cv.lin$var1.pred ~ cv.lin$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_lin <- lm(cv.lin$var1.pred ~ cv.lin$observed)
abline(lm_lin, col = "red", lwd = 2)

r2_lin = summary(lm_lin)$r.squared
rmse_lin = hydroGOF::rmse(cv.lin$var1.pred, cv.lin$observed)

## 2) Modelo exponencial
cv.exp <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = exp.fit)
summary(cv.exp)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred           var1.var         observed        residual       
 Min.   : 0.06051   Min.   : 12.38   Min.   : 0.00   Min.   :-15.0221  
 1st Qu.: 3.53749   1st Qu.: 18.23   1st Qu.: 2.20   1st Qu.: -2.5773  
 Median : 5.69650   Median : 24.71   Median : 5.00   Median : -0.5739  
 Mean   : 7.01900   Mean   : 31.82   Mean   : 7.34   Mean   :  0.3210  
 3rd Qu.: 9.60323   3rd Qu.: 38.94   3rd Qu.:10.30   3rd Qu.:  2.6853  
 Max.   :20.22214   Max.   :107.38   Max.   :36.00   Max.   : 30.4535  
     zscore              fold      
 Min.   :-2.50740   Min.   : 1.00  
 1st Qu.:-0.46053   1st Qu.: 8.25  
 Median :-0.11106   Median :15.50  
 Mean   : 0.02631   Mean   :15.50  
 3rd Qu.: 0.45814   3rd Qu.:22.75  
 Max.   : 4.65631   Max.   :30.00

plot(cv.exp$var1.pred ~ cv.exp$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_exp <- lm(cv.exp$var1.pred ~ cv.exp$observed)
abline(lm_exp, col = "red", lwd = 2)

r2_exp = summary(lm_exp)$r.squared
rmse_exp = hydroGOF::rmse(cv.exp$var1.pred, cv.exp$observed)

## 3) Modelo Gaussiano
cv.gau <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = gau.fit)
summary(cv.gau)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred         var1.var         observed        residual        
 Min.   : 0.471   Min.   : 12.83   Min.   : 0.00   Min.   :-15.95270  
 1st Qu.: 2.869   1st Qu.: 14.63   1st Qu.: 2.20   1st Qu.: -2.33169  
 Median : 6.791   Median : 19.70   Median : 5.00   Median : -0.46101  
 Mean   : 7.297   Mean   : 29.96   Mean   : 7.34   Mean   :  0.04298  
 3rd Qu.: 9.971   3rd Qu.: 40.45   3rd Qu.:10.30   3rd Qu.:  1.88845  
 Max.   :21.153   Max.   :105.30   Max.   :36.00   Max.   : 30.93771  
     zscore               fold      
 Min.   :-2.765848   Min.   : 1.00  
 1st Qu.:-0.464753   1st Qu.: 8.25  
 Median :-0.101036   Median :15.50  
 Mean   : 0.004199   Mean   :15.50  
 3rd Qu.: 0.248281   3rd Qu.:22.75  
 Max.   : 4.660612   Max.   :30.00

plot(cv.gau$var1.pred ~ cv.gau$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_gau <- lm(cv.gau$var1.pred ~ cv.gau$observed)
abline(lm_gau, col = "red", lwd = 2)

r2_gau = summary(lm_gau)$r.squared
rmse_gau = hydroGOF::rmse(cv.gau$var1.pred, cv.gau$observed)

## 4) Modelo Wave
cv.wav <- krige.cv(acumulado_24h ~ x + y, locations = pluvio, model = wav.fit)
summary(cv.wav)

Object of class SpatialPointsDataFrame
Coordinates:
        min     max
x  634918.7  687966
y 7450222.8 7475508
Is projected: NA 
proj4string : [NA]
Number of points: 30
Data attributes:
   var1.pred         var1.var        observed        residual       
 Min.   : 1.979   Min.   :14.15   Min.   : 0.00   Min.   :-11.4223  
 1st Qu.: 4.483   1st Qu.:20.05   1st Qu.: 2.20   1st Qu.: -4.7160  
 Median : 6.881   Median :27.94   Median : 5.00   Median : -0.5548  
 Mean   : 7.024   Mean   :28.33   Mean   : 7.34   Mean   :  0.3161  
 3rd Qu.: 9.391   3rd Qu.:35.76   3rd Qu.:10.30   3rd Qu.:  2.8856  
 Max.   :13.217   Max.   :48.97   Max.   :36.00   Max.   : 31.7935  
     zscore              fold      
 Min.   :-1.94021   Min.   : 1.00  
 1st Qu.:-0.93925   1st Qu.: 8.25  
 Median :-0.13349   Median :15.50  
 Mean   : 0.03151   Mean   :15.50  
 3rd Qu.: 0.52837   3rd Qu.:22.75  
 Max.   : 5.09689   Max.   :30.00

plot(cv.wav$var1.pred ~ cv.wav$observed, cex = 1.2, lwd = 2)
abline(0, 1, col = "lightgrey", lwd = 2)
lm_wav <- lm(cv.wav$var1.pred ~ cv.wav$observed)
abline(lm_wav, col = "red", lwd = 2)

r2_wav = summary(lm_wav)$r.squared
rmse_wav = hydroGOF::rmse(cv.wav$var1.pred, cv.wav$observed)

# Criando uma tabela das estatística de R2 e RMSE
df.r2 <- data.frame(r2_lin, r2_exp, r2_gau, r2_wav)
df.rmse <- data.frame(rmse_lin, rmse_exp, rmse_gau, rmse_wav)
tabela <- data.frame(cbind(t(df.r2), t(df.rmse)))
colnames(tabela) <- c("R2", "RMSE")
rnames <- gsub("r2_", "", rownames(tabela))  # remove o prefixo r2 dos nomes das linhas
rownames(tabela) <- rnames  # substitui o nome das linhas simplificadas na tab original
tabela

              R2     RMSE
lin 0.1482468633 6.882779
exp 0.0765441167 7.496573
gau 0.0680566292 7.723425
wav 0.0001772582 7.991183

Criando os grids do contorno da cidade do Rio de Janeiro para intermpolação

# Importando o contorno do Rio
contorno.rio <- shapefile("dados/chuva_rio/MUNIC_2K_2022_IPP_SIRGAS.shp")

# Criando grade para interpolacao com a resolucao de 50m
r <- raster(contorno.rio, res = 50)

# Criando um objeto formato raster
rp <- rasterize(contorno.rio, r, 0)

# Trasnsformando o objeto raster no formato SpatialPixelsDataFrame
grid <- as(rp, "SpatialPixelsDataFrame")
sp::plot(grid)

Krigagem

# Colocando os dados de chuva e o grid na mesma projecao
sp::proj4string(pluvio) = CRS(proj4string(contorno.rio))

mapa_chuva_lin <- krige(acumulado_24h ~ 1, pluvio, grid, model = lin.fit)

[using ordinary kriging]

mapa_chuva_exp <- krige(acumulado_24h ~ 1, pluvio, grid, model = exp.fit)

[using ordinary kriging]

mapa_chuva_gau <- krige(acumulado_24h ~ 1, pluvio, grid, model = gau.fit)

[using ordinary kriging]

mapa_chuva_wav <- krige(acumulado_24h ~ 1, pluvio, grid, model = wav.fit)

[using ordinary kriging]

plot.lin <- plot(mapa_chuva_lin, main = "Modelo Linear")

plot.exp <- plot(mapa_chuva_exp, main = "Modelo Exponencial")

plot.gau <- plot(mapa_chuva_gau, main = "Modelo Gaussiano")

plot.wav <- plot(mapa_chuva_wav, main = "Modelo Wave")

Auto Krige

# Modelando
auto.krige = autoKrige(acumulado_24h ~ x + y, pluvio, grid, model = "Exp")

[using universal kriging]

summary(auto.krige)

krige_output:
Object of class SpatialPixelsDataFrame
Coordinates:
        min       max
x  623533.4  695333.4
y 7446574.5 7483374.5
Is projected: TRUE 
proj4string :
[+proj=utm +zone=23 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs]
Number of points: 481750
Grid attributes:
   cellcentre.offset cellsize cells.dim
s1          623558.4       50      1436
s2         7446599.5       50       736
Data attributes:
   var1.pred           var1.var           var1.stdev     
 Min.   :-0.01083   Min.   :  0.07308   Min.   : 0.2703  
 1st Qu.: 1.50803   1st Qu.: 14.16177   1st Qu.: 3.7632  
 Median : 5.25560   Median : 22.76226   Median : 4.7710  
 Mean   : 6.66402   Mean   : 24.90904   Mean   : 4.7804  
 3rd Qu.:11.20533   3rd Qu.: 31.22818   3rd Qu.: 5.5882  
 Max.   :35.96625   Max.   :148.38326   Max.   :12.1813  

exp_var:
  np      dist    gamma dir.hor dir.ver   id
1  8  2254.389 11.19978       0       0 var1
2 27  4384.026 13.98185       0       0 var1
3 35  6073.421 36.23544       0       0 var1
4 54  8869.824 53.69680       0       0 var1
5 50 12016.434 41.83955       0       0 var1
6 47 15019.442 52.23264       0       0 var1
7 52 18659.164 85.79906       0       0 var1

var_model:
  model    psill    range
1   Nug   0.0000     0.00
2   Exp 176.6964 33077.08
Sums of squares betw. var. model and sample var.[1] 0.0003263469

# Validação cruzada
auto.krige.cv <- autoKrige.cv(acumulado_24h ~ x + y, pluvio, model = "Exp")
summary(auto.krige.cv)

            [,1]   
mean_error  0.36   
me_mean     0.04904
MAE         4.544  
MSE         57.79  
MSNE        1.789  
cor_obspred 0.2754 
cor_predres -0.3772
RMSE        7.602  
RMSE_sd     1.022  
URMSE       7.594  
iqr         4.636

Convertendo para o formato raster - auto krige

raster_chuva <- raster(auto.krige$krige_output)
plot(raster_chuva)

Caso queira salvar a imagem raster em um arquivo formato geotiff para ler em algum SIG por exemplo:

# Exportando o objeto da imagem
writeRaster(raster_chuva, filename = "dados/chuva_rio/chuva_auto.tiff", format = "GTiff",
    overwrite = T)

# Importando de volta para o R
raster_chuva <- raster("dados/chuva_rio/chuva_auto.tiff")

Fazendo o mapa interativo com as estações

# Importando os dados referente as estatções pluviométricas
estacoes.sf <- read_sf("dados/chuva_rio/Estac_C3_B5es_Alerta_Rio.shp")

# Convertendo UTM para Lat Long das estacoes
estacoes.longlat <- st_transform(estacoes.sf, "+proj=longlat +ellps=WGS84 +datum=WGS84")
estacoes.longlat$coords <- st_coordinates(estacoes.longlat)
estacoes.longlat$X <- estacoes.longlat$coords[, 1]
estacoes.longlat$Y <- estacoes.longlat$coords[, 2]

# Importando a malha de bairros
bairros.sf <- read_sf("dados/chuva_rio/BAIRROS_2K_2022_IPP_SIRGAS.shp")

# Convertendo UTM para Lat Long a malha dos bairros
bairros.longlat <- st_transform(bairros.sf, "+proj=longlat +ellps=WGS84 +datum=WGS84")

# Convertendo o raster da chuva para lat long
raster_chuva_longlat <- projectRaster(raster_chuva, crs = CRS("+proj=longlat +datum=WGS84"))

# Definindo Paleta de cores da superfície interpolada da chuva
pal <- colorNumeric(c("#000066", "#00c8f8", "#F0E68C", "#FFFF00", "#FF8C00"), values(raster_chuva_longlat),
    na.color = "transparent", reverse = T)

# Construindo o mapa interativo via leaflet

leaflet(data = estacoes.longlat, options = leafletOptions(attributionControl = FALSE)) |>
    # addTiles() |>
addProviderTiles("CartoDB.Positron", group = "Ruas") |>
    addProviderTiles("Esri.WorldImagery", options = providerTileOptions(opacity = 0.7),
        group = "Satélite") |>
    addProviderTiles(providers$CartoDB.Voyager, group = "Voyager") |>
    addProviderTiles(providers$Stamen.Toner, group = "Toner") |>
    setView(lng = -43.42, lat = -22.9, zoom = 10.4) |>
    # addProviderTiles(providers$CartoDB.Voyager) |>
addMarkers(~X, ~Y, popup = ~as.character(est), label = ~as.character(est), group = "Estações") |>
    ############## Polígonos dos Bairros ################
addPolygons(data = bairros.longlat, weight = 3, color = "darkblue", smoothFactor = 1,
    fill = FALSE, labelOptions = labelOptions(style = list(`font-weight` = "normal",
        padding = "3px 8px"), textsize = "13px", direction = "auto"), group = "Bairros") |>
    ########## Adicionando o raster #########################
addRasterImage(raster_chuva_longlat, colors = pal, opacity = 0.8, group = "Chuva: 1 semana") %>%
    leaflet::addLegend(pal = pal, values = values(raster_chuva_longlat), title = "Chuva Acumulada - 1 semana",
        group = "Chuva: 1 semana") |>
    ############## Controle das layers (botoes) ################
addLayersControl(baseGroups = c("Voyager", "Ruas", "Satélite", "Toner"), overlayGroups = c("Estações",
    "Bairros", "Chuva: 1 semana"), options = layersControlOptions(collapsed = FALSE),
    position = "bottomleft") |>
    ########## Desabilitando os grupos ################
hideGroup(group = c("Bairros"))

Bibliografia sugerida

Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds). Análise Espacial de Dados Geográficos. Brasília, EMBRAPA, 2004.
Interactive Spatial Data Analysis by Trevor C. Bailey , Anthony C. Gatrell Routledge, 1995
Applied Spatial Statistics for Public Health Data; Lance A. Waller, Carol A. Gotway Wiley-Interscience 1St ed. 2004
Applied Spatial Data Analysis with R; Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer; Edição: 2nd ed. 2013
Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow. Geocomputation with R, 2021.
Spatial Statistics Workbook of Department of Criminology at the University of Manchester by Reka Solymosi and Juanjo Medina Crime Mapping in R
Spatial Data Science with R Spatial Data Science with R
GeoComputation and Spatial Analysis practicals by Lex Comber On-line Book
GEOG5917 Big Data & Consumer Analytics - RStudio Practicals by Lex Comber On-line Book

Especialização em Estatística Aplicada (PPGEAD): Módulo de Estatística Espacial

Introdução à Análise Estatística Espacial

O que é Análise Estatística Espacial ?

Origem da Estatística Espacial

Objetivos da Estatística Espacial

Dependência Espacial ou Autocorrelação Espacial

Tipologia dos Dados Espaciais

Padrões Pontuais

Exemplos de aplicação

Estimativas de Kernel (ou Mapas de Calor)

Como funciona:

Geoestatística

Estrutura dos dados geoestatísticos

Exemplos de aplicação

Semivariograma: a base da análise espacial

Dados de Área

Exemplos de aplicação

Mapa Temático

Matriz de Vizinhança

Moran Global, Moran Local e Lisa Map

🧮 Moran Global

🧮 Moran Local (LISA)

🗺️ LISA Map

Modelos de Regressão Espacial

i) Modelos Globais:

ii) Modelos Locais:

Aplicação I: Utilizando a biblioteca tmap para construção rápida de mapas temáticos

Aplicação II: Baixando e construindo mapas a partir da biblioteca geobr

Bibliotecas que serão utilizadas:

Primeiro, vamos fazer um gráfico apenas com as geometrias.

Aplicação III: Dengue em Dourados/MS - Parte 1: Análise exploratória

Biliotecas do R que serão utilizadas

Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS

Fazendo um join com as tabelas com os setores censitários + população

Lendo e plotando os casos de dengue georreferenciados em Dourados/MS

Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS

Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.

Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)

Representando os casos e o lixo de forma interativa com o a biblioteca tmap

Convertendo o dado de pontos (padrão pontual) para dados de área

Plotando o mapa temático dos casos por setor censitário

Plotando o mapa temático dos pontos de coleta de lixo por setor censitário

Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário

Plotando a distribuição da incidência em Dourados/MS

Kernel por atributos

Extraindo os centróides dos polígonos em Dourados/MS

Colocando os pontos no formato sp

Construindo o kernel por atributo da taxa de detecção

Construindo a matriz de vizinhança para verificar a autocorrelação espacial

Verificando a malha de conectividade da vizinhança de Dourados/MS

Obtendo a correlação da taxa de incidência de dengue Dourados/MS

Plotando o correlograma

Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.

Moran Local (Lisa Map) da taxa de incidência Dourados/MS

LISA Map

Aplicação IV: Dengue em Dourados/MS - Parte 2: Modelagem (Modelos Linear, CAR e GWR)

Ajustando o modelo de regressão linear simples.

Checando os residuos para verificar a presença de autocorrelação.

Ajustando o modelo CAR (Spatial Error Model)

Checando os residuos para verificar a presença de autocorrelação

Ajustando o modelo GWR (Geographically Weighted Regression)

Verificando a distribuição dos coeficientes de regressão para a variável lixo

Verificando a distribuição dos localR2

Incorporando alguns parâmetros de saída do modelo na tabela setor.sf

Mapa dos coeficientes de regressão para a variável lixo

Checando os residuos para verificar a presença de autocorrelação para o modelo GWR.

Mapeando os coeficientes de regressão para a variável lixo por significancia através do teste de wald

Mapa dos coeficientes de determinação regionalizados (\(R^2\) local).

Verificando a distribuição dos preditos.

Mapeando os valores observados e preditos dos modelos ajustados

Comparando a distribuição dos resíduos dos modelos ajustados

Mapeando os resíduos dos modelos ajustados

Aplicação V: Geoestatística com dados de pluviosidade na cidade do Rio de Janeiro/RJ

Biliotecas do R que serão utilizadas

Importando a tabela com a chuva acumulada média de 7 dias das últimas 24hs e 96hs das estações pluviométricas da cidde do Rio de Janeiro

Análise Gráfica Descritiva

Transformar os dados em um objeto espacial do R

Análise Gráfica Descritiva com os dados espaciais

Modelando o variograma experimental (ou empírico)

Ajustando o semivariograma teórico

Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando `st_intersection`.

Representando os casos e o lixo de forma interativa com o a biblioteca `tmap`